rm(list = ls())

# Load required libraries
library(tidyverse)
library(haven)

# Set working directory 
setwd("~/Dropbox/Whitten_Kagalwala/TS/JOP/Acceptance/KW_replication/WLL Bounds Approach")

#---- Coverage data preparation ----
bounds <- read_dta("bounds_coverage.dta")

# Load lrm type1 data from "UnitRoot_Independent" folder
lrm <- read_dta("unitroot_independent.dta")
lrm$cov_adl22 <- 1-lrm$lrm_adl22_t1
lrm$cov_adl33 <- 1-lrm$lrm_adl33_t1
lrm <- select(lrm, c(t, cov_adl22, cov_adl33))

# Merge
data <- left_join(bounds, lrm, by = "t")
data <- select(data, c(t, ftr_adl, ftr_adl22, ftr_adl33, cov_adl22, cov_adl33))

# Reshape
data <- gather(data, Method, coverage, ftr_adl:cov_adl33)
data$Method[data$Method == "ftr_adl"] <- "Bounds-ADL(1,1)"
data$Method[data$Method == "ftr_adl22"] <- "Bounds-ADL(2,2)"
data$Method[data$Method == "ftr_adl33"] <- "Bounds-ADL(3,3)"
data$Method[data$Method == "cov_adl22"] <- "General-ADL(2,2)"
data$Method[data$Method == "cov_adl33"] <- "General-ADL(3,3)"
data$Method <- as.factor(data$Method)
data$t <- as.factor(data$t)

#---- Figure 8a, Bounds Independently Generated Unit Roots Coverage ----
ggplot(data, aes(t, coverage, fill = Method)) +
  geom_hline(yintercept = 0.95, linetype= "dashed", color=gray(1/2)) + 
  geom_bar(position = "dodge", stat = "identity") + 
  labs(x = "T", y = "", fill = "", title = "Coverage of LRM") + ylim(0,1) +
  theme_bw() +
  theme(plot.title = element_text(size = 10), legend.position = "bottom", legend.box = "horizontal",
                   legend.text = element_text(size = 6), axis.title = element_text(size = 10)) +
  guides(fill = guide_legend(nrow = 2, byrow = TRUE)) +
  scale_fill_grey()
ggsave("coverage_bounds.pdf", width = 20, height = 20, units = "cm")
ggsave("coverage_bounds.eps", device = "eps", width = 20, height = 20, units = "cm")

#---- Power data preparation ----
bounds_adl <- read_dta("bounds_ADL.dta")

# Load lrm type1 data
lrm_adl <- read_dta("adl.dta")
lrm_adl <- lrm_adl %>%
  filter(phi_y == 0.6 & phi_x == 0.6) %>%
  dplyr::select(t, lrm_adl22_power, lrm_adl33_power)

# Merge
data_adl <- left_join(bounds_adl, lrm_adl, by = "t")
data_adl <- dplyr::select(data_adl, c(t, rej_adl, rej_adl22, rej_adl33, lrm_adl22_power, lrm_adl33_power))

# Reshape
data_adl <- gather(data_adl, Method, power, rej_adl:lrm_adl33_power)
data_adl$Method[data_adl$Method == "rej_adl"] <- "Bounds-ADL(1,1)"
data_adl$Method[data_adl$Method == "rej_adl22"] <- "Bounds-ADL(2,2)"
data_adl$Method[data_adl$Method == "rej_adl33"] <- "Bounds-ADL(3,3)"
data_adl$Method[data_adl$Method == "lrm_adl22_power"] <- "General-ADL(2,2)"
data_adl$Method[data_adl$Method == "lrm_adl33_power"] <- "General-ADL(3,3)"
data_adl$Method <- as.factor(data_adl$Method)
data_adl$t <- as.factor(data_adl$t)

#---- Figure 8b, Bounds ADL DGP Power ----
ggplot(data_adl, aes(t, power, fill = Method)) +
  geom_hline(yintercept = 1, linetype= "dashed", color=gray(1/2)) + 
  geom_bar(position = "dodge", stat = "identity") + 
  labs(x = "T", y = "", fill = "", title = "Power of LRM") + ylim(0,1) +
  theme_bw() +
  theme(plot.title = element_text(size = 10), legend.position = "bottom", legend.box = "horizontal",
        legend.text = element_text(size = 6), axis.title = element_text(size = 10)) +
  guides(fill = guide_legend(nrow = 2, byrow = TRUE)) +
  scale_fill_grey()
ggsave("power_bounds.pdf", width = 20, height = 20, units = "cm")
ggsave("power_bounds.eps", device = "eps", width = 20, height = 20, units = "cm")

#---- Misspecification data preparation ----
bounds_PA <- read_dta("bounds_ADL.dta")
bounds_PA <- bounds_PA %>%
  dplyr::select(t, lrm_pa, ftr_pa, inc_pa, rej_pa)

# Reshape
bounds_PA <- gather(bounds_PA, Rate, value, ftr_pa:rej_pa)
bounds_PA$Rate[bounds_PA$Rate == "ftr_pa"] <- "Fail to Reject"
bounds_PA$Rate[bounds_PA$Rate == "inc_pa"] <- "Inconclusive"
bounds_PA$Rate[bounds_PA$Rate == "rej_pa"] <- "Reject"
bounds_PA$Rate <- as.factor(bounds_PA$Rate)
bounds_PA$t <- as.factor(bounds_PA$t)

#---- Figure 8c, Bounds ADL DGP Misspecification PA  ----

ggplot(bounds_PA, aes(t, value, fill = Rate)) +
  geom_bar(position="stack", stat="identity") +
  labs(x = "T", y = "Rate", fill = "", title = "Model Misspecifcation and WLL Bounds",
       caption = "Partial adjustment model estimated on ADL(1, 1) DGP") + ylim(0,1) +
  theme_bw() +
  theme(plot.title = element_text(size = 10), legend.position = "bottom", legend.box = "horizontal",
        legend.text = element_text(size = 6), axis.title = element_text(size = 10)) +
  guides(fill = guide_legend(nrow = 1, byrow = TRUE)) +
  scale_fill_grey() + 
  theme(plot.caption.position = "plot",
        plot.caption = element_text(hjust = 0))
ggsave("misspec_bounds.pdf", width = 20, height = 20, units = "cm")
ggsave("misspec_bounds.eps", device = "eps", width = 20, height = 20, units = "cm")

